Convergence of all-order many-body methods: coupled-cluster study for Li 
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We present and analyze results of the relativistic coupled-cluster calculation of energies, hyperfine 
constants, and dipole matrix elements for the 2s, 2pi/ 2 , and 2p 3 / 2 states of Li atom. The calculations 
are complete through the fourth order of many-body perturbation theory for energies and through 
the fifth order for matrix elements and subsume certain chains of diagrams in all orders. A nearly 
complete many-body calculation allows us to draw conclusions on the convergence pattern of the 
coupled-cluster method. Our analysis suggests that the high-order many-body contributions to 
energies and matrix elements scale proportionally and provides a quantitative ground for semi- 
empirical fits of ab inito matrix elements to experimental energies. 

PACS numbers: 31.15.ac, 31. 15. am, 32.10.Fn, 32.70.Cs 

Many-body perturbation theory (MBPT) is a ubiqui- 
tous tool in atomic and nuclear physics and quantum 
chemistry. Yet its order-by-order convergence has been 
found to fail in several systems (see, e.g., [l|). To circum- 
vent this drawback, one usually employs all-order meth- 
ods which implicitly sum most important classes (chains) 
of diagrams to all orders of MBPT. Even in this case, as 
we illustrate here with a nearly-complete solution of the 
many-body problem for Li, the saturation with respect to 
a systematic addition of the all-order chains may reveal 
a non-monotonic convergence. In other words, including 
increasingly complex (and computationally more expen- 
sive) chains does not necessarily translate into a better 
accuracy. 

While such a convergence pattern may seem discour- 
aging, we find that the high-order many-body contribu- 
tions to energies and matrix elements vary proportionally 
as the all-order formalism is augmented with increasingly 
complex chains of diagrams. We explain this dependence 
by the similarity of self-energy contributions to both en- 
ergies and matrix elements and provide a quantitative 
ground for semi-empirical fits. This is especially valuable 
for atomic systems, where high-accuracy experimental 
data for energies are available, while the matrix elements 
have a relatively poor accuracy. In some cases, e.g., in 
parity violation, the matrix elements of the weak interac- 
tion are not known experimentally at all, while they need 
to be computed to a high precision [2, 3]. Although the 
semi-empirical fits have been used before [3J, [j] , the va- 
lidity of such scaling has not been rigorously established. 
Here, based on a nearly complete many-body calculation, 
we are able to address this question. 

We solve the many-body problem for the three-electron 
Li atom. Here the availability of both high-accuracy 
variational Hylleraas and experimental data makes the 
analysis of minute high-order MBPT effects plausible. 
Our calculations are complete through the fourth order 
of MBPT for energies and through the fifth order for ma- 
trix elements. Additionally certain classes of diagrams 
are summed to all orders using the coupled cluster (CC) 
method. The previous CC-type formulations for Li [fj|, Q 



were complete only through the second order for energies 
and the third order for matrix elements. 

We consider Li as a univalent atom and choose the 
lowest-order Hamiltonian to include the relativistic ki- 
netic energy operator of electrons and their interac- 
tions with the nucleus and the V N ~ l Dirac-Hartree-Fock 
(DHF) potential. The single-particle orbitals and ener- 
gies £i are found from the set of the frozen-core DHF 
equations. Using the DHF basis, the Hamiltonian reads 
(up to an energy offset) 

H = H + G = J2 SiN[a\ai] + - ^ gwN[a\a]aia k ] . 

i ijkl 

(i) 

Here Hq is the one-electron lowest-order Hamiltonian, G 
is the residual Coulomb interaction, a] and a, are the cre- 
ation and annihilation operators, and N[- • ■ ] is the nor- 
mal product of operators with respect to the core quasi- 
vacuum state |0 C ). Indices i,j, k and I range over all pos- 
sible single-particle orbitals, and gijki are the Coulomb 
matrix elements. 

We are interested in obtaining the exact many-body 
state \*$/ v ) that is seeded from the lowest-order DHF state 

l*i 0) )=4|0c): 
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(2) 



where $7 (yet to be found) is the so-called wave opera- 
tor Q. In the CC method the MBPT diagrams are re- 
summed to all orders and one introduces the exponential 
ansatz for the wave operator 



n = N[exp(K)) = 1 + K + ^N[K 2 } 



(3) 



where the cluster operator K is expressed in terms of con- 
nected diagrams of the wave operator f2. The operator 
K is decomposed into cluster operators (K) n combining 
n simultaneous excitations of core and valence electrons 
from the reference state W v ) to all orders of MBPT. For 



the three-electron Li, the exact cluster operator reads 
K = S c + D c + S v + D v +T V = (4) 



the double-headed arrow representing the valence state. 
S v and D v (S C ,D C ) are the valence (core) singles and 
doubles, and T v are the valence triples. This exhausts 
the entire excitation basis for the three-electron Li (e.g., 
there are no core triples). 

In the previous CC work for Li [5[ the expansion ([J| 
was truncated at the S and D excitations and the CC 
equations contained only terms linear in the CC ampli- 
tudes (SD method). The full CC study involving S and 
D excitations was carried out in [8( ; we will refer to it as 
the CCSD approximation. Our present treatment is nat- 
urally labeled as CCSDvT to emphasize our additional 
inclusion of the valence triple excitations. 



cUtiUdiLJ e 



S V [T V ] D V [T V ] SE V [T V ] Z[T V ] 






T„ \D C 



m 



T V [D V 



T„ \T„] 





T v [D c (g> D c ] 



T v [S v <g> D c ] 



T„ \D„ ® D r . 



FIG. 1: Representative diagrams for various classes of con- 
tributions of the valence triples. The wavy lines represent 
the residual Coulomb interaction and those capped with the 
heavy square represent a one-body interaction (e.g., hyperfine 
interaction) . 



A set of coupled equations for the cluster operators 
(K) n ( (K c ) 1 = S c , (K v ) 1 = S v , etc.) may be found 
from the Bloch equation [7f specialized for univalent sys- 
tems n 

(s v -H )(K c ) n = {QGn} conncctcdn , 
(s v + 5E V H ) (K v ) n = {QG fi} connected ,„ , (5) 

where the valence correlation energy 

6E V = (*W\GSl\*W) (6) 

and Q = 1 — \^ v )($? v | is the projection operator. 



Below we present a topological structure of the CC 
equations for the cluster amplitudes in the CCSDvT ap- 
proximation. The resulting equations for the core cluster 
amplitudes S c and D c are the same as in the CCSD ap- 
proximation [10| | and we do not repeat them here. Repre- 
sentative diagrams involving triples are shown in Fig. [T] 
The structure of the valence singles equation is 



[H , S v ] + 6E V S V = CCSD + S V [T V ] 



(7) 



where notation like ^[T^] stands for the effect of the va- 
lence triples (T v ) on the r.h.s. of the equation for valence 
singles (S v ). Here [Hq,S v ] is a commutator, and SE V is 
the valence correlation energy, 



6E V = SEccsd + SE V [T V ] . 



(8) 



where SEqcsd is obtained within the CCSD approach 
and SE V [T V ] is due to the valence triples. The equation 
for the valence doubles reads 

-[H , D v ] + 5E V D V « CCSD + D V [T V ] . 

Here we discarded contribution D v [S c ® T v ] which stands 
for a nonlinear contribution resulting from a product of 
clusters S c and T v . For the valence triples we obtain 



-[H ,T V ]+SE V T V « T V [D C 
T V [D C ®D V ] + T V [D C 



■T V [D V ]+T V [T V ] + 
D c ) + T v [S v ® D c 



with the discarded terms of higher order in G. Our ap- 
proximation subsumes the entire set of fourth order dia- 
grams for SE V . This is a substantial improvement over 
both the SD and the CCSD method which are complete 
only through the second order of MBPT. 

Solution of the CCSDvT equations provides us with 
the wavefunctions and the correlation energies. With 
the obtained wavefunctions we compute the matrix el- 
ements. The relevant CCSDvT formalism is presented 
in Refs. [H], El- -"- n addition to the well-explored SD 
contributions [5|, our formalism includes contributions 
from the valence triples and also "dressing" of matrix 
elements. The dressing arises from re-summing non- 
linear contributions to the atomic wavefunctions @ 
in expressions for matrix elements and, in particular, 
guarantees that the important chain of random-phase- 
approximation diagrams is fully recovered in all orders 
of MBPT [ll(. In addition we incorporated all contri- 
butions that are quadratic in valence triples. Overall the 
calculations of matrix elements are complete through the 
fifth order of MBPT and incorporate certain classes of di- 
agrams summed to all orders. 

Our numerical calculations are based on our previ- 
ous CCSDvT code [10], with the addition of the entire 
set of the non-linear CCSD contributions documented in 
Ref. Q. The important new additions are the effects 
of triples on triples T„[T„] and the leading-order non- 
linear terms on the r.h.s. of the triples equations. We also 
employed efficient dual-kinetic-balance basis sets, as de- 
scribed in Ref. [12| • The basis set included partial waves 



£ = — 6 for the S and D amplitudes and £ = — 4 for the 
T v amplitudes. Our final results include extrapolation for 
an infinitely large basis. Since the hyperfine interaction 
occurs at small distances, due to the uncertainty relation, 
one has to keep orbitals with high excitation energies in 
the basis. 



TABLE I: Contributions to removal energies of 2s, 2p 1/ / 2 , and 
2p3/2 states for Li in cm -1 in various approximations. 





2si/ 2 


2pi/ 2 


2P3/2 


DHF 


43087.3 


28232.9 


28232.3 


ASD 


405.8 


352.0 


351.9 


ACCSD 


-6.8 


-5.0 


-5.0 


A{T V [D V + D c ]) a 


2.8 


2.6 


2.6 


A(T V [T V }) 


2.8 


3.3 


3.3 


A(T v [NL]) a 


-1.1 


-1.1 


-1.1 


corrections 6 


-3.3(5) 


-0.7(5) 


-0.4(5) 


Total 


43487.5 


28584.1 


28583.7 


Experiment [13] 


43487.2 


28583.5 


28583.2 




Other CC works 




SD+MBPT-III [6] 


43487.5 


28581.9 


28581.5 


CCSD [14] 


43483 


28567 





T v [D v + D c ] = T v [D v ] + T v [D c 

T v [NL] = T v [D c <g> D v ] + T v [D c <g> D c ] + T v [S v <g> D c ] 

'includes basis set, recoil, Breit, and QED corrections. Error 

bar is due to basis extrapolation 




MBPT complexity — ► 

FIG. 2: (Color online) Convergence pattern of the CCSDvT 
method as a function of MBPT complexity for the HFS con- 
stant of the ground state of Li. Breit, QED, recoil and basis 
set corrections are included in all theoretical values. 
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TABLE II: Contributions to the magnetic-dipole hyperfine 
structure constants A of 2s, 2p!/ 2 , and 2p 3 / 2 states for 7 Li 
(I = 3/2, fi — 3.256427(2)) in MHz in various approxima- 
tions. 



CCSD 



-1.5 - § 





2Si/ 2 


2P1/2 


2P3/2 


DHF 


284.35 


32.295 


6.457 


ASD 


117.68 


13.622 


-9.474 


ACCSD 


-1.79 


-0.233 


0.176 


A dressing 


-0.40 


-0.039 


0.031 


A(T V [D V + D C }) 


1.25 


0.218 


-0.175 


A(T„[T„]) 


0.28 


0.058 


-0.031 


A(T„[NL]) 


-0.04 


-0.010 


0.002 


corrections™ 


0.33(3) 


0.046(6) 


-0.026(1) 


Total 


401.66 


45.958 


-3.041 


Experiment 


401. 75. . b 


45.914(25)" 


-3.055(14)" 


Experiment 




46.010(25) d 





"includes basis set, recoil, Breit, and QED corrections. Error 

bar is due to the basis extrapolation. 

"401.7520433(5) Schlecht and McColm [l]|; c Orth et al. []J| 

d Walls et al. El; 



In Tables fl] and |TT] we present calculated energies and 
magnetic-dipole hyperfine structure (HFS) constants A. 
In these tables the entries are ordered by increasing 
MBPT complexity of the calculations. A denotes a dif- 
ference from the preceding entry due to extra classes of 
the diagrams included at that level of approximation. 
For example, the entry ACCSD is obtained by taking 
a difference between the CCSD and the SD results. We 
include Breit, QED, and recoil corrections in our final 



FIG. 3: (Color online) Variations from the the final CCS- 
DvT values of A2s and E23 in different approximations, e.g., 
AXsd = Xsd—XccsdvT- The variations in matrix elements 
and energies are correlated and exhibit linear dependence. 



result. For energies they were adopted from Ref. [6J, for 
A2S from pjl, and for A2 P from |19fl . The basis set ex- 
trapolation was carried out in the conventional manner 
(see, e.g., |6j). The error bar is estimated as a half of the 
basis-set extrapolation correction. The HFS constants 
were computed using finite nucleus and uniform magne- 
tization. 

The correlation contributions follow a similar pattern 
in all these cases. We illustrate the convergence of the 
all-order method in the case of the HFS constant for 
the ground state in Fig. [5] Here the experimental un- 
certainty is about 1 ppb so that the deviation from the 
experimental value is an indication of the theoretical ac- 
curacy. The dominant correlations are recovered at the 
SD level, which captures all third-order diagrams and re- 
sults in a 0.2% theoretical accuracy. The inclusion of the 
CCSD non-linear effects and the dressing leads to a worse 
agreement with the experiment (-0.4%). An inclusion of 
the leading valence triples returns the agreement to the 
0.1% level. The addition of the higher-order T V [T V ] effect 



results in an almost perfect agreement with the exper- 
iment. Finally, T„[NL] diagrams provide only a minor 
correction. The final result is complete through the fifth 
order and agrees with the experiment at the 0.02% level. 

Finally, we present the computed reduced ma- 
trix elements of the electric-dipole operator. We 
obtain in the CCSDvT approximation with the 
dressing {2si/ 2 \\D\\2pi/ 2 ) = 3.31633(7) |e|ao and 
(2s 1/2 \\D\\2 P3/2 ) = 4.6901(1) |e|a . These values also in- 
clude basis extrapolation and are complete through the 
fifth order of MBPT. The error bars here correspond to a 
half of the basis extrapolation correction. Again, the con- 
vergence with respect to the addition of higher-order dia- 
grams follows the same pattern as for the HFS constants 
(Fig. [2]) and energies. To facilitate a comparison with 
the previous high-accuracy variational studies we form 
the oscillator strength / for the 2s — 2p transition. Our 
result, / = 0.74686, is smaller than the non-relativistic 
variational value (201 by 0.01%. The size and the sign of 
the difference is consistent with the expected difference 
due to relativistic effects. 

The accuracy and the completeness of our calculations 
allow us to make the following observations. 

Accuracy. The CCSDvT method improves the accu- 
racy over the previous less complete CC-type calcula- 
tions. For example, for energies the overall agreement 
stands at a few 0.1 cm -1 while the SD method is accu- 
rate to a few cm -1 . Similarly, there is an order of magni- 
tude improvement in the accuracy of computing A 2s over 
that of the SD approximation. The remaining differences 
with the experiment are likely due to higher-order dia- 
grams discarded in our scheme; these are consistent with 
the size of the T„[NL] effect. 

Convergence. In the absence of general theorems on 
convergence of MBPT, the present work provides an em- 
pirical proof that the CC method converges for Li. For 



all the computed properties the saturation of the method 
with respect to adding increasingly complex classes of di- 
agrams is not monotonic, as illustrated in Fig. [2 The 
empirical conclusion for other, more complex, univalent 
atoms is that both the non-linear CCSD effects and the 
valence triples have to be treated simultaneously [2| ■ 

Correlation between corrections to the energy and to 
the matrix elements. There is a strong link between 
the convergence patterns for energies and for matrix el- 
ements. This is illustrated in Fig. [3J the deviations of 
the A 2s and E 2s from the final CCSDvT values follow 
roughly a linear law. The data for matrix elements does 
not include dressing. A similar pattern is observed for 
other matrix elements as well. Such a linear dependence 
is due to the effect of self-energy (Brueckner) correction. 
This dominant chain of diagrams is presented in both 
matrix elements and energies. For example, for triple ex- 
citations, the corrections 5V, [X 1 ,,] and SE V [T V ] arise from 
the same diagram and the modification of singles due to 
triples propagates into the calculation of the matrix ele- 
ment. Similar scaling ideas were used earlier [3, |4j to fit 
low-order results to higher orders, but never rigorously 
tested. Of course, as apparent from Fig. [3J the linear 
scaling is only approximate and can be used in the semi- 
empirical fits only to a certain accuracy. For example, the 
self-energy corrections do not affect "dressing" of matrix 
elements which contributes at a sizable 0.1% level to the 
A 2s constant. Neither can it capture the distinctively- 
different QED corrections to the energies and matrix el- 
ements. 
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